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ABSTRACT 

Supermassive black hole binaries (SMBHBs) should be an inevitable consequence of the hierarchical 
growth of massive galaxies through mergers and the strongest sirens of gravitational waves (GWs) in 
the cosmos. Yet, their direct detection has remained elusive due to the compact (sub-parsec) orbital 
separations of gravitationally bound SMBHBs. Here, we exploit a theoretically predicted signature 
of an SMBHB in the time domain: periodic variability caused by a mass accretion rate that is 
modulated by the binary’s orbital motion. We report our first significant periodically varying quasar 
detection from the systematic search in the Pan-STARRSl (PSI) Medium Deep Survey. Our SMBHB 
candidate, PSO J334.2028-I-0I.4075, is a luminous radio-loud quasar at z = 2.060, with extended 
baseline photometry from the Gatalina Real-Time Transient Survey, as well as archival spectroscopy 
from the EIRST Bright Quasar Survey. The observed period (542 ±15 days) and estimated black 
hole mass (log(MBH/Al0) = 9.97 ± 0.50), correspond to an orbital separation of Schwarzschild 
radii (~ O.OOGtp ooa Pc), assuming the rest-frame period of the quasar variability traces the orbital 
period of the binary. This SMBHB candidate, discovered at the peak redshift for SMBH mergers, 
is in a physically stable configuration for a circumbinary accretion disk and within the regime of 
GW-driven orbital decay. Our search with PSl is a benchmark study for the exciting capabilities 
of LSST, which will have orders of magnitude larger survey power and will potentially pinpoint the 
locations of thousands of SMBHBs in the variable night sky. 

Subject headings: quasars: general — surveys 


1. INTRODUCTION 

The expectation for the existence of supermassive 
black hole binaries (SMBHBs) in galaxy nuclei is sup¬ 
ported by two well-established properties of galax¬ 
ies: (I) high spatial resolution observations of nearby 
galaxies have demonstrated that SM BHs are ubiqui¬ 
tous in the centers of galaxy bulges (|Magorrian et al.l 
[TM^ with masses tightly correla ted with the mass 
and s t ructure of their host g a laxies (iFerra rese fc MerrittI 
l200rt iGebhardt et al.l l200(l iGraham et al.11200IH . and 
(2) galaxies in a AGDM Universe build up their struc¬ 
ture hierarchically through mergers (e.g. iSoringel et al.l 
1200511 . When two galaxies merge, their SMBHs will sink 
to the center through dynamical friction and through 
three-body interactions with stars and viscous exchange 
of angular momentum with circumbinary gas form a 
gravitationally bound binary that eventually coalesces 
due to the radiatio n of gravitational waves (GWs; 
iBegelman et al.iri980H . 

Recent progress has been made in the detection of 
“dual active galactic nuclei (AGNs),” double active nu¬ 
clei in assumed merged galaxy system s with kiloparsec- 
scale separations (|Komossa et alJl^003l : IComerford et al.l 
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1200911 . These dual AGNs, while a product of a galaxy 
merger, are not yet gravitationally bound, and thus are 
not necessarily fated to coalesce. A true SMBHB be¬ 
comes gravitationally bound on the scale of parsecs, 
which, beyond our Local Group of galaxies, is well below 
the angular resolving power of the most powerful cur¬ 
rent, or even future, telescopes. However, several promis¬ 
ing candidates have been identified indirectly via spec¬ 
troscopy: quasars with offset and/or drifting broad-line 
peaks attributed to a broad-lin e region in orbit around 
an SMBH’s bina r y companion ([Boroson fc Laued 120091 
IDotti et all 120091 fearrows et al.il20IIIl . However, alter¬ 
native scenarios have been proposed that do not require 
an SMBHB, in cluding double-peaked lines from a single 
accretion disk (iChornock et al.l[2010ll . 

A promising observational signature of SMBHBs is 
their variable accretion luminosity. One of the first sub¬ 
parsec SMBHB cand idates, OJ287, was ident ified by its 
variability behavior (ILehto fc Valton^ 1199611 . OJ287 is 
a quasar that demonstrates regular optical outbursts on 
a timescale of 12 yr that has been modeled as the result 
of a secondary SMBH companio n passing through th e 
primary SMBH’s accretion disk (IValtonen et aLlfeoPSH . 
Such a configuration should be rare since the secondary 
BH’s orbital axis must be highly misaligned with the 
primary BH’s accretion disk axis in order for it to pass 
through its disk. A more generic signature of an SMBHB 
is likely to be related to accretion through its circumbi¬ 
nary disk. 
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In a gas-rich galaxy merger, strong gravitational 
torques drive gas i nward, triggering bot h star formation 
and BH accretion (|HoDkins et al.l 1200611 . In particular, 


hydrodynamical simulations of circumbinary disks show 
that accretion via “hot streams” onto the BHs is strongly 
modulated by the bi nary’s orbital motion for mass ra- 


tios of 0.05 < q < 1 (iMacFadven & Milosavlievid 

2008 

Shi et al.l 120121: INoble et al.l 120121: iD’Orazio et al.f 

2013 

Gold et al.l I2014D. Simulations (ID’Orazio et al.l 

2013 

Farris et al.ll2015D also detect a t 6 Urs timescale orig- 


inating from a surface density “lump” just outside the 
central cavity of the circumbinary disk, which is a per¬ 
sistent but secularly evolving feature. 

While working on our paper, we became aware of the 
report of PG 1302-102, a periodically variable quasar 
discove red by the Catalina Real-Time Transient Survey 
(CRTS: [Graham et al.ll201^ . It is a 15 magnitude quasar 
at z = 0.2784, varying at the 0.14 mag level with a pe¬ 
riod of 5.2 ± 0.2 yr, with good sampling over 1.8 cy¬ 
cles and extended archival data going back 20 yr. Their 
physical interpretation for its variability is an SMBHB 
{\og{M/M q) ~ 8.5, a ^ 0.01 pc), its luminosity being 
modulated due to either a precessing jet or an overden¬ 
sity (“hot spot”) in the inner edge of its circumbinary 
disk. 

In this Letter, we present our most significant detec¬ 
tion from the systematic search for periodically varying 
quasars in the Pan-STARRSl Medium Deep Survey (PSl 
MDS) field MD09. PSO J334.2028-k01.4075 is a radio- 
loud quasar at z = 2.060 with archival spectroscopy from 
FQBS and extended baseline photometry from CRTS. 
The 8.5 yr baseline of the PSU-CRTS light curve is well 
described by a simple sinusoid, consistent with theoret¬ 
ical simulations for the modulated accretion rate in a 
0.05 < q < 0.25 mass-ratio SMBHB. We use the rest- 
frame period and virial black hole mass estimate to infer 
an orbital separation of the binary that is in the GW- 
driven regime. 

2. THEORETICAL PREDICTIONS 

The dynamics of a gravitationally bound SMBHB sys¬ 
tem can be described by Kepler’s third law: 

{whi) (lo^) ■ 

where i?s is the Schwarzschild radius, Rg = a is 

the separation betw een the BHs; and M is the total 
mass of the system. iHaiman et al.l (120091) calculate the 
minimum survey area to detect a statistically significant 
sample of quasars powered by SMBHBs as a function 
of variable magnitude depth, assuming reasonable val¬ 
ues for the quasar luminosity function, quasar lifetime 
{tg = 10^ yr), the Eddington fraction (/Edd = 0.3), and 
the fractional variability amplitude {^f/f = 0.1). The 
variability detection threshold we have achieved in MD09 
((j4|) corresponds to a A/// > 0.1 sensitivity for point 
sources brighter than m ~ 21 mag, and t hus a variable 
magni tude of 23.5 mag. At this depth, IHaiman et al.l 
(|2009D require an area of ~100 deg^ to yield a sample 
of over 100 SMBHBs; an excellent match to the area of 
PSl MDS (80 deg^). Furthermore, the baseline (4.2 yr) 
and cadence (3 days) of PSl MDS make us sensitive to 


timescales for which BHs (with M > 10^ Mq) are in the 
GW-driven regime of orbital decay. 

3. PAN-STARRSl MEDIUM DEEP SURVEY 

The Panoramic Survey Telescope and Rapid Response 
System (Pan-STARRS) is a wide-held imaging system 
designed for dedicated survey observations on a 1.8m 
telescope on Haleakala, Hawaii, w ith a 1.4 Gigapixe l 
camera and a 7 deg^ held of view (|Kaiser et al.l 1201^ . 
The PSl telescope is operated by the Institute for As¬ 
tronomy (IfA) at the University of Hawaii and has just 
completed over 4 yr of operation in 2014 March. We 
present data from the Medium Deep Survey (MDS), a 
deep, multi-epoch survey of 10 circular helds distributed 
across the sky, each ^ 8 deg^ in size, whose daily observ¬ 
ing cadence in hve hlters is excellent for studying persis¬ 
tent variable sources, including quasars. The PSl MDS 
cadence of observation cycles through the gpi, rpi, ipi, 
and Zpi bands every three nights, with observations in 
the ypi band close to the full Moon. Due to the poorer 
time sampling of the ypi observations, we do not use 
them in this analysis. 

4. ENSEMBLE PHOTOMETRY 

We began our systematic search for SMBHB candi¬ 
dates among color-selected quasars in the PSl MD09 
held. This is the hrst MD held that was made available 
to the PSl Science Gollaboration in the Pan-STARRS 
Science Interface (PSl) online database. In order to 
maximize our sensitivity to intrinsic variability, we hrst 
applied the technique of differential ensemble photom¬ 
etry. This technique is able to correct for local sys¬ 
tematic errors due to variable atmospheric conditions 
by co mparing a target object with nearby non-variable 
stars (iHonevcuttl 119921: iBhatti et al.ll20inl i. We created 
a color-selected reference star sample and quasar sample 
by cross-matching point sources (m < 23 mag) in the 
MD09 held with a custom catalog extracted from full- 
survey deep stacks from PSl MDS in the gpi, rpi, ipi, 
zpi, and ypi bands, as well as from observations with 
the Canada-France-Hawaii Telescope (GFHT) in the u 
band (S. Heinis et al. 2015, in preparation). We con- 
verted all magnitu des to the SDSS photometric system 
(iTonrv et al.ll201^ to take advantage of the SDSS color 
selection of stars and quasars already available in the lit¬ 
erature (| Schmidt et al.l 120101 : iSesar et fo] l2007t) . Figure 
[T] shows the color-color diagrams of the point sources in 
MD09 selected as quasars and non-variable stars. This 
query resulted in 8158 reference stars and 316 quasars, 
each with an average of 350 detections in four hlters. 

We m odihed the en s emble photometry software devel¬ 
oped by IBhatti et al.l (|20inD for SDSS to the PSl data 
format and ran it on the reference stars. In Figure [1] 
(bottom right panel), we plot the “corrected” magnitude 
error as a function of mean magnitude compared with the 
“raw” values before ensemble photometry. The ensem¬ 
ble photometry reduces the measured errors signihcantly, 
lowering the error hoor from 0.045 to 0.025 mag, and re¬ 
sulting in a 2tT variability threshold of 0.05 mag on the 
bright end to 0.34 mag on the faint end. 

5. SELECTING PERIODIC QUASAR CANDIDATES 

We then applied ensemble photometry to the 316 color- 
selected quasars and hagged quasars as variable based 
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on their magnitude error relative to their neighbors of a 
similar brightness; we set 2 (t as our criterion for vari¬ 
ability and required a variability flag in at least two 
filters. This selection yielded 168 variable quasars in 
MD09. Among these color-selected variable quasars, 
we searched for potential periodic signatures using the 
Lomb-Scargle (LS) periodogram, a Fourier analysis tech¬ 
nique of unevenly spaced data with no ise (iLombllTQTGl : 
iScargld 119821 : iHorne fc Baliunad I1986I1 . For Nq data 
points in the time series spanning a total length of T in 
units of MJD, we sampled the periodogram at the num- 
ber of recommend e d inde pendent frequencies (VQ, from 
iHorne fc Baliunai (|1986ll . from 1/T to No/{2T) (which 
would be the Nyquist frequency if data were evenly sam¬ 
pled), resulting in a frequency resolution in the peri¬ 
odogram of A/ = {No/2 — l)/{TNi) ~ 2 X 10“^ d“^. 
When identifying periodic sources, we took advantage of 
the redundancy of PSl MDS monitoring in four filters 
(ffpi ?’pi ipi zpi), each with a slightly different observ¬ 
ing cadence due to weather and scheduling constraints 
to help filter out false detections from aliasing by requir¬ 
ing that periodogram peaks be coherent across multiple 
filters; 40 of the 168 variable quasars survived this test. 

6. PERIODIC QUASAR CANDIDATE PSO J334.2028 
-1-01.4075 

Among the candidate periodic quasars from our pe¬ 
riodogram analysis, here we focus on our most signif¬ 
icant detection, PSO J334.2028-b01.4075 (J2000). In 
Figure [H we show its periodogram in four PSl filters, 
with the strongest peak marked with a dashed line. 
We fold each filter light curve on this period (Figure 
ED, and measure the scatter of the residuals from the 
best-fit sine curve {cr)- The error of the periodogram 
peak frequency {6f) and the signal-to-noise ratio of the 
peak power (^) can then be calculated from tr^ and the 
amplitude of the signal Aq (IHorne &: Baliunaa 1 198611 as 
6f = (which gives us an error on the detected 

period of SP = Sf/p) and ^ = Ag/{2a'/.), respectively. 

The resulting average period across all four filters is 
P = 541.8±15.3 days, with the highest signal-to-noise ra¬ 
tio in the gpi filter with ^ = 3.19, and periodogram peaks 
in all four filters well above a 1.5x 10“^^ false-alarm prob¬ 
ability (corresponding to lOcr) plotted with a dotted line 
in Figured The PSl data cover 2.6 cycles, just shy of the 
“rule of thumb” number of cycles (t hree) for a p eriodic 
variation to be apparent to the eye l|Pressl [19781) . From 
our Monte Carl o simulations of 1 000 stochastic Damped 
Random Walk ([Kelly et al.ll2009ll light curves, we find a 
false periodic detection rate of 6.3% using our selection 
criteria from §5. We further disfavor a false alarm from 
stochastic quasar variability since the 0.6% of the sim¬ 
ulations that successfully mimic the periodic timescale 
of our candidate have short-timescale variances a factor 
of > 2 larger than expected for the quasar’s luminosity 
and inferred black hole mass (§7). Note that there is a 
secondary peak in the gpi and rpi periodograms that, 
if real and not an artifact from the PSl data sampling, 
is at twice the primary peak frequency, a signature of 
0.05 < q < 0.25 mass-ratio SMBHBs, which show an 
accretion rate m odulation most close ly described by a 
simple sinusoid (iD’Orazio et al.l 1201311 . The amplitude 
of PSO J334.2028-1-01.4075’s sinusoidal modulation in¬ 


creases with decreasing wavelength, consistent with the 
exponential dependence on wa velength found in previous 
quasar variabilitY stu dies fe.g. iVanden Berk et al.ll2004 
iMacLeod et al.ll2010l l. 

PSO J334.2028-1-01.4075 is a radio-loud quasar (FBQS 
J221648.7-1-012427) with an ar chival spectrum fro m the 
FIRST Bright Quasar Survey (iBecker et al.ll200lll . We 
are also fortunate that this ca ndidate has an arc hival V- 
band light curve from CRTS (iDrake et al.ll2009ll . which 
we use to test the persistence of the periodic variation 
over an extended baseline of 8.5 yr (corresponding to 5.7 
cycles). To compare to the CRTS light curve, we con- 
yert the PSl qpi-b and light curves to the SDSS system 
([Tonrv et al.l 1201^ and then to the Johnson V magni¬ 
tude using_the_photometric transformation for quasars 
from lJester et al.l (|2005h and an average gpi-rpi = -bO.lO 
mag. We had to apply an additional offset of —0.17 mag 
to the pseudo-K PSl magnitudes in order to match the 
average of the CRTS K-band data. Though the pho¬ 
tometric errors are relatively large, the CRTS measure¬ 
ments are consistent with those of PSl during their over¬ 
lap (Figure HI), and have residuals over the entire CRTS 
baseline from the PSl-fitted sinusoidal model that are 
Gaussian with a a — 0.17 mag that is comparable to the 
mean photometric error of 0.18 mag. 

7. PHYSICAL INTERPRETATION 

We use the width of the quasar’s C IV line and its 
nearby continuum lumin osity to make a virial estimate o f 
the black hole mass from iVestergaard fc Peterson! (1200611 : 


/Mbh\ 




/FWHM(C IV) 
lOOOkm/s 


Y( ) 
J V lO'^'^erg/s/ 


, 0.53i 


+6. 


( 2 ) 

where A = 1350 A . The FBQS spectrum, though not 
publicly available in an electronic format, was mea¬ 
sured with a ruler to determine FA,obs(I350A (U-z)) 
~ 8.5 X 10“^^ ergs s“^cm“^ A“^, and FWHM (CIV 
AI550) ~ 200 A. The CIV line is symmetric in shape, 
and its width corresponds to a velocity in the rest-frame 
of ~ 12,650 km s“^. We co rrect for a Galactic extinctio n 
of E{B — V) = 0.0406 mag (ISchlaflv &: Finkb eineill20II[l . 
using the extinction law from iCardelli et al.l (|1989ll , to 
find LA.em(1350l) = FA,obsl0^^/2 ®(l -b z) - 9.5 x 
10^^ erg s“^A“^, where is the luminosity distance 
assuming iJo = 70 km/s/Mpc, VLm = 0.3, Da = 0.7, 
and k = 0, and A\ = 0.1790. This results in a black 


hole mass of 


(^) 


= 9.97, with a scatter from the 


uncertainty in the relation of 0.5dex. Applying a mean 


quasar bolometric cor rection at I350A of BC=3.81 from 
iRichards et al.l ([200611 . one gets a bolometric luminosity 
of Lboi = XL\BC = 4.9 X 10"^® erg s“^. Note that this 
object is also radio loud, with a radio luminosity at the 
rest- frame frequen c y of 5 GHz of log(L/i(erg s“^)) =32.8 
from IBecker et aP (200111 . 

Assuming the rest-frame period Prest = Pobs/(l + z) 
is on the order of the orbital timescale of the SMBHB, 
with a caveat that in addition to a strong dependence 
on mass ratio, there are a range of theoretical predic¬ 
tions for translating P-Rst. to torb (e.g. iNoble et al.ll2012L 
ID’Orazio et al.ll2013ll . we then calculate the orbital sep- 
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aration of the binary to be 7'!l^Rs O.OObto^ooa pc), 
securely placing it in the gravitationally bound regime 
of a physically viable SMBHB system — a circumbinary 
accretion disk system capable of maintaining a central 
cavity, stable to gravitational fragmentati on, and in the 
regime of orbital decay driven by GWs (iHaiman et alJ 

120091 : li^csis et al.ll2012l : lD’Orazio et al.ll2013ir ' Also note 

that since the viscous time scales as one could ex¬ 
pect to be able to detect modulations in the accretion 
rate fed by the streams in the circumbinary disk cav¬ 
ity, without being washed out by viscous processes i n 
the “mini disks” around each BH (|Roedig et alJl201^ . 
Remarkably, the rest-frame inspiral time for the binary 
is tinsp = = 7.0yr(a/7R.)4(M/109-97Me)-3 

for q = 0.25, where /r = opening up the possi¬ 

bility for detecting the decay of the orbital period (P) 
with future monitoring, as well as providing a promising 
target for direct GW detection for pulsar timing arrays 
(|Sesana et al.l[200^ . 


8. DISCUSSION AND CONCLUSIONS 


We present the most statistically significant periodi¬ 
cally variable quasar candidate from our search in PSl 
MD09, PSO J334.2028-1-01.4075, a radio-loud quasar at 
z = 2.060. We combine an estimate of its black hole mass 
with its variability timescale (assuming Prest iorb) 
to find orbital parameters consistent with model pre¬ 
dictions of a stable accreting SMBHB s ystem with a 
0.05 < q < 0.25 in the GW-driven regime ( Haiman et al.l 
[ 200 ^ . 


The redshift of this SMBHB candidat e coincides with 
the peak epoch for SMBHB mergers (IVolonteri et al.l 
l2003f) . and its large mass {M « IO^^Mq) makes it favor¬ 
able for detection in the GW-driven regime, due to the 
strong dependence on M of the res idence time at a give n 
orbital separation in units of Rs (IHaiman et al.l 1200911 . 
Like the CRTS SMBHB candidate PG 1302-102 reported 


bv I Graham et al.l (1201511 . our SMBHB candidate is also 
a radio-loud quasar. However, given the shorter rest- 
frame period of our candidate of 0.5 yr (versus 4 yr in 
PG 1302-102), it is even more unlikely for its variabil¬ 
ity to be driven b y jet precession, e ither originating from 
a single SMBH (iLu fc Z hod l2005f) or a binary SMBH 
(|Lobanov fc Rolandll2005ll . 

This pilot program in PSl MD09 is a promising start 
to our systematic search for periodic variability signa¬ 
tures of SMBHBs amongst the expected ~ 1000 variable 
quasars in the full ^ 80 deg^ PSl MDS. At the start 
of the next decade (~ 2023), the Ly ge Synoptic Survey 
Telescope (LSST) (llvezic et al.ll2008ll will probe a volume 
several thousand times larger than PSl MDS, yielding 
tens of millions of quasars, and potentially thousands of 
SMBHBs periodically varying on the timescale of years, 
fated to coalesce. 
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Fig. 1. — Color-color diagrams used to select point sources in MD09 for the 8158 point sources in the reference star sample (red) and 
316 point sources in the quasar sample (blue). Photometry is measured from the CFHT+PSl catalog and converted to the SDSS system. 
Dashed lines show the color selection boundaries. The stellar color-color selection box was chosen to avoid RR Lyrae stars, which are 
intrinsically variable. Bottom right panel: observed standard deviation of the reference sample of non-variable stars before and after 
applying the technique of ensemble photometry (open circles and filled squares, respectively), compared to the Poisson error expected from 
the reported PSl flux errors (black dashed lines). 
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Fig. 2.— Our automated LS periodogram routine selects periodically variable candidates by requiring that the strongest peak is detected 
at the same frequency in at least three filters. This quasar candidate, PSO J334.2028+01.4075, was selected through this method and had 
the periodogram peak with the highest signal-to-noise ratio of all of our candidates. The dashed lines mark the strongest peak in each 
filter. The dotted line corresponds to a false-alarm probability of 1.5 x 10“^^, or 10 a. 
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Fig. 3.— Upper panels: sinusoidal fit to the folded PSl light curve of PSO J334.2028+01.4075 in four filters, with the error bar indicating 
the typical photometric error for an object of similar brightness in that filter. The period corresponding to the peak of the periodogram and 
its error bar, the amplitude of the fitted sine wave, and the signal-to-noise ratio of the peak power are each labeled. Lower panel: sinusoidal 
fit plotted over the complete PSl light curves in the gpi rpi ipi and zpi bands. The data used to create this figure are available. 
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Fig. 4. — PSl light curve (asterisks) converted to the V band to compare to the archival CRTS light curve (dots with error bars). The 
CRTS data points are binned in one-day intervals, with the error bars measured from the standard deviation in the bin and not including 
data points with a photometric error greater than 0.25 mag or nights with less than three measurements. This results in 34/113 nights of 
data being thrown out. The CRTS measurements are overall consistent with the PSl light curve and the sine fit to the PSl light curve 
(dashed curve). 



